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We present a systematic approach for calculating higher-order derivatives of smooth functions 
on a uniform grid using Pade approximants. We illustrate our findings by deriving higher-order 
approximations using traditional second-order finite-differences formulas as our starting point. We 
employ these schemes to study the stability and dynamical properties of K{2, 2) Rosenau-Hyman 
(RH) compactons including the collision of two compactons and resultant shock formation. Our 
approach uses a differencing scheme involving only nearest and next-to-nearest neighbors on a 
uniform spatial grid. The partial differential equation for the compactons involves first, second 
and third partial derivatives in the spatial coordinate and we concentrate on four different fourth- 
order methods which differ in the possibility of increasing the degree of accuracy (or not) of one 
of the spatial derivatives to sixth order. A method designed to reduce roundoff errors was found 
to be the most accurate approximation in stability studies of single solitary waves, even though all 
derivates are accurate only to fourth order. Simulating compacton scattering requires the addition 
of fourth derivatives related to artificial viscosity. For those problems the different choices lead to 
different amounts of "spurious" radiation and we compare the virtues of the different choices. 
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I. INTRODUCTION 

Since their discovery by Rosenau and Hyman in 
1993 [1], compactons have found diverse applications in 
physics in the analysis of patterns on liquid surfaces [ ] , 
in approximations for thin viscous films [ ], ocean dy- 
namics [4], magma dynamics [ , ], and medicine [ ]. 
Compactons are also the object of study in brane cos- 
mology [ ] as well as mathematical physics [ , ], 
and the dynamics of nonlinear lattices [11-14] to model 
the dispersive coupling of a chain of oscillators [^^-1"]. 
Multidimensional RH compactons have been discussed 
in [18, 19]. Recently, compact structures also have been 
studied in the context of a Klein-Gordon model [ , ]. 
A recent review of nonlinear evolution equations with co- 
sine/sine compacton solutions can be found in Ref. 

Compactons represent a class of traveling-wave solu- 
tions with compact support resulting from the balance of 
both nonlinearity and nonlinear dispersion. Compactons 
were discovered by Rosenau and Hyman (RH) in the pro- 
cess of studying the role played by nonlinear dispersion in 
pattern formation in liquid drops using a family of fully 
nonlinear Korteweg-de Vries (KdV) equations [ ] , 



where u = u{x^ t) is the wave amplitude, x is the spatial 
coordinate and t is time. 

RH called these solitary waves compactons, and 
Eq. (1.1) is known as the K(l^p) compacton equation. 
The RH compactons have the remarkable soliton prop- 
erty that after colliding with other compactons they 
reemerge with the same coherent shape. However, unlike 
the soliton collisions in an integrable system, the point 
where the compactons collide is marked by the creation 
of low-amplitude compacton-anticompacton pairs [ ] . 

The RH generalization of the KdV equation (1.1) is 
only derivable from a Lagrangian in the K{1^ 1) case. 
Hence, in general, Eq. (1.1) does not exhibit the usual 
energy conservation law. Therefore, Cooper, Shepard 
and Sodano [ ] proposed a different generalization of 
the KdV equation based on the first-order Lagrangian 



Lir,s) = J 



r{r — 1) 



+a{^xyi^^^f dx. (1.2) 



Ut + (U% + K)a 



0, 



(1.1) 
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We note that the set (/,p) in Eq. (1.1) corresponds to the 
set (r — l,s + l) in Eq. (1.2). Since then, various other La- 
grangian generalizations of the KdV equation have been 
considered [9zL-9R]. With the exception of Ref. ['^•^], the 
structural stability of the resulting compactons was stud- 
ied solely using analytical techniques such as linear sta- 
bility analysis [ ], and an exhaustive numerical study 
of the stability and dynamical properties of these com- 
pacton solitary waves is needed. 

In general, the numerical analysis of compactons is a 
difficult numerical problem because compactons have at 
most a finite number of continuous derivatives at their 
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edges. Unlike the comp actons derived from the La- 
grangian (1.2), the RH compact ons have been the object 
of intense numerical study using pseudospectral meth- 
ods [ , ], finite-element methods based on cubic B- 
splines [ , ] and on piecewise polynomials discontinu- 
ous at the finite element interfaces [31], finite-differences 
methods [zz, ^z-^o], methods of lines with adaptive mesh 
refinement [36, 37], and particle methods based on the 
dispersive- velocity method [38]. 

Both the pseudospectral and finite-differences meth- 
ods require artificial dissipation (hyperviscosity) to sim- 
ulate interacting compactons without appreciable spuri- 
ous radiation. The RH pseudospectral methods use a 
discrete Fourier transform and incorporate the hyper- 
viscosity using high-pass filters based on second spatial 
derivatives. Using this approach RH showed successfully 
that compactons collide without any apparent radiation. 
However, because the pseudospectral methods explicitly 
damp the high-frequency modes in order to alleviate the 
negative effects due to the high-frequency dispersive er- 
rors introduced by the lack of smoothness at the edges 
of the compacton, the pseudospectral approach is not 
suitable for the study of high-frequency phenomena, and 
the usability of filters themselves has been called in the 
question [ ]. In turn, finite-differences methods usually 
incorporate the artificial dissipation via a fourth spa- 
tial derivative term. However, in the absence of high- 
frequency filtering, these methods are marred by the ap- 
pearance of spurious radiation [33] . This radiation prop- 
agates both backward and forward and has an amplitude 
smaller by a few orders of magnitude than the compacton 
amplitude. Its numerical origin can be identified by a grid 
refinement technique. 

Recently, Rus and Villatoro [22, 33-35] introduced 
a discretization procedure for uniform spatial grids 
based on a Fade approximant-like [ ] improvement of 
finite-differences methods. As special cases, this ap- 
proach can be used to obtain the familiar second-order 
finite-differences methods and the fourth-order Fetrov- 
Galerkin finite-element method developed by Sanz-Serna 
and co-authors [29, 30]. Given the involved character 
of the Fetrov-Galerkin approach based on linear inter- 
polants described in Ref. [30], we believe the work by 
Rus and Villatoro (RV) lends itself to further scrutiny. 

In this paper we present a systematic derivation of the 
Fade approximants [39] intended to calculate derivatives 
of smooth functions on a uniform grid by deriving higher- 
order approximations using traditional finite-differences 
formulas. Our derivation recovers as special cases the 
Fade approximants first introduced by Rus and Villa- 
toro [22, 33-35]. We illustrate our approach for the par- 
ticular case when second-order finite-differences formulas 
are used as the starting point to derive at least fourth- 
order accurate approximations of the first three deriva- 
tives of a smooth function. This approach is equivalent 
to deriving the best differencing schemes involving only 
nearest and next-to-nearest neighbors on a uniform grid. 
We apply these approximation schemes to the study of 



stability and dynamical properties of K{p^p) Rosenau- 
Hyman compactons. This study is intended to estab- 
lish the baseline for future studies of the stability and 
dynamical properties of L(r, s) compactons, which fea- 
ture higher-order nonlinearities and terms with mixed- 
derivatives that are not present in the K{p^p) equations. 
Hence, the numerical analysis of the properties of the 
L(r, s) compactons of Eq. (1.2) is expected to be consid- 
erably more difficult. 

This paper is outlined as follows. In Sec. II, we show 
that the approximation schemes discussed by RV [33, 34] 
can be identified as special cases of a systematic improve- 
ment scheme that uses Fade approximants to derive at 
least fourth-order accurate approximations for the first 
three spatial derivatives Um , Um , and Um (we have in- 
troduced the spatial discretization Xm = mh^ u{x) -^ Um 
and the roman numeral superscript denotes the order of 
spatial derivative at x^) by starting with second-order 
finite-differences approximations. In general, one can 
begin with finite-differences approximations of any arbi- 
trary even order, and improve upon these by at least two 
orders of accuracy by using suitable Fade approximants. 
In Sec. Ill we discuss several special cases: Three of these 
cases describe approximation schemes that mix fourth- 
order accurate approximations for two of the derivatives 
Um ^ Uni\ and Um with a sixth-order accurate approx- 
imation for the third one. We also discuss the case of 
the "optimal" fourth-order approximation scheme. In the 
latter, all three derivatives are fourth-order accurate, but 
all the coefficients entering the Fade approximants have 
values that result in a reduction of decimal roundoff er- 
rors. In Sec. IV we apply the above four approximation 
schemes to study the stability and dynamical properties 
of i^(2, 2) Rosenau-Hyman compactons. We conclude by 
summarizing our main results in Sec. V. 



II. FADE APPROXIMANTS 

To begin, we consider a smooth function u{x)^ defined 
on the interval x G [0,L], and discretized on a uniform 
grid, Xm = mh^ with m = 0, 1, • • • , M, and h = L/M. 
Fade approximants of order k of the derivatives of u{x) 
are defined as rational approximations of the form 



u^ =4^u^ + 0{Ax'' 



(2.1) 



-^ =||iy«™ + ^(^^')' (2.2) 



(Hi) 



^u^ + 0{Ax'^) 



HE) 



w^i") = 3S w™ + ©(Ax'^) 



HE) 
where we have introduced the shift operator, E^ as 



(2.3) 
(2.4) 



E^Um = Um+k- (2.5) 



In this language, the second-oidei accurate approxima- 
tion of derivatives based on finite-differences correspond 
to the Pade approximants given by [ ] 



AiiE) 
B,{E) 
Ci{E) 
Vi{E) 



1 



2Aa; 
1 

1 

2Aa;3 
1 



AxH- 



E-E-'^ 

E-2 + E-^ 

E^ -2E + 2E-^ - E-^ 
E"^ -4E + 6- AE-^ + ^-2 



(2.6) 

(2.7) 

(2.8) 

, (2.9) 



and Ti{E) = 1. We note that even- and oc^c^-order 
derivatives require approximants that are symmetric and 
antisymmetric in E^ respectively. 

We also note that although all four operators, Ai{E)^ 
Bi{E), Ci{E), and Vi{E), lead to second-order accu- 
rate numerical approximations, the derivatives ^Xm and 
Um involve the subset of grid points {xm-,Xm±i-,Xm±2}^ 
whereas the derivatives Um and vr^ involve only the 
subset of grid points {xm^Xm±i}- Therefore, it is possi- 
ble to design a numerical scheme that improves the order 
of approximation of the derivatives Um and Um by in- 
corporating the additional grid points, {^^^±2} (see Ap- 
pendix A). 

It is more challenging, however, to find a consis- 
tent approach that improves the order of approxima- 
tion of all four lowest-order derivatives without extend- 
ing the set of grid points. We will show next that the 
Pade-approximant approach described here, allows us to 
provide a consistent approach involving only the grid 
points {xm^ Xm±i^Xm±2} that includes three of these four 
derivatives. 

In the following we will use extensively the Taylor ex- 
pansion of u{x) around x^, i.e. 

Um+k = u{xm + kAx) = Um ^ u^^ {kAx) (2.10) 



Xii) 



k^Ax^ 



(Hi) 



JfAa 
6~ 



(iv) 



k^Ax^ 
24 



(^) k^Ax^ (^^) k^Ax^ (^^^) k^Ax^ 
+ ^m ^20 ^ "^ 720 ^ "^ 5040 ^ 

The following two relationships follow immediately: 



(iS'^ 



H/ I Urn — '^m-\-k ~r Ufn—k — ^ ^m 



(2.11) 



u^::h'Ax' + v^^ ^^^ ^ .A-) ^^ 



12 



360 



and 



(^^ - E-^) Um = Um+k - Um-k = 2u^^kAx (2.12) 



E ■ J u 

(iii) '^ ^^ 



/.(-) 



k^Ax^ 
60 



(vii) 



k^Ax^ 
2520 



To obtain a fourth-order accurate approximation of 
the derivatives, we can either begin by improving the 



third-order derivative, u^ -, or the fourth-order deriva- 
tive, Um • Unfortunately, we cannot improve both these 
derivatives at the same time. Because in the compacton- 
dynamics problem [1, 29-34], the fourth-order derivative 
enters only through the artificial viscosity term needed 
to handle shocks, we chose to improve the approxima- 
tion corresponding to the third-order derivative, Um • 

In the following we derive the operators A2{E)^ B2{E)^ 
C2{E)^ and V2{E), corresponding to the new fourth-order 
accurate Pade approximants. 



A. Third-order derivatives 



Using Eqs. (2.3) and (2.12), we obtain 

Ci{E)Um 



(2.13) 



— n,iiii) 



,iv) 



Ax^ 



(vii) 



.^Ax' 



40 



(ix) 



17Ax^ 
12096 



or 



,(**«) — 



Ci{E)u, 



t^^^O{Ax'). 



(2.14) 



To eliminate the dependence on Ax^, we consider a linear 
combination of the second-order approximations of i^m 
on the same subset of grid points, {xm^Xm±i,Xm±2}- 
This can be achieved by introducing an operator, J^{E), 
symmetric in E, such that 



TiE) u, 



{Hi) 



1 r 



(£2 



such that 

T{E)ul^"^=C^{E)ur, 
Using Eq. (2.11), we obtain 



b{E + E-'^)+c 



0{Ax'' 



"to ' 

(2.15) 
(2.16) 



2u(l") + «(;) Aa;2 + «("") ^ + ■ • 



4Aa;'* 



■ cul 



' r 

(2.17) 



Requiring that this approximation is fourth order or bet- 
ter, we obtain the system of equations 



a-2b-c 
a- 4h 



2, 
16. 



and its solution can be parameterized as: 

a = 4r, 6 = r-4, c = 2(r + 3) 
It follows that we can write 
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(2.19) 

(2.20) 
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24 
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Hence, we have 
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1 1\ A 
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43 
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V60 rJ 

1\ Ax^ 



x^ 



(2.22) 



24 



(9(Ax«) . 



For r integer and r > 5, we obtain solutions with a, 6, 
and c positive integers. 



and we can write 



with 

A2{E) = 



A2{E) 

T{E) "^"^ ^^ V30 r. 



^ V30 r; 



1\ Ax^ 



Vl05 4r/ 6 



4 

(9(Ax^) 



24Ax 



£;^ + lOE - lOE-^ - E-^ 



(2.31) 



(2.32) 



First-order derivatives 



Second-order derivatives 



Next, we calculate the corresponding Fade approxi- 
mants for the first-order derivative, Um- We consider 



.(0 



A{E) 
HE) 



0{Ax^) , 



(2.23) 



with J-{E) given by (2.21), and require that the order of 
the approximation is fourth order or better. Therefore, 
A{E) must be an operator antisymmetric in E. 
By definition, we introduce 



A{E)u, 



1 



and solve 



aAx 



{E 



2 ^-2^ 



(E-E-'] 



T{E) uW = A{E) Um + 0{/^x^) . 



(2.24) 



(2.25) 



We have 

A{E) Ur. 
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8Aa;2 
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™ 315 
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60 
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(2.26) 



To satisfy the requirement of a fourth-order accurate ap- 
proximation for the first-order derivative Um\ we solve 
the system of equations 



a - 2/3 = 4 , 
3a - 4/3 = 32 , 

and obtain the solution 

a = 24 , /3 = 10 . 

This gives 

,Ax2 



(2.27) 
(2.28) 



(2.29) 



A2{E)Ur, 



= ./(O 



u, 



{Hi) ' 



(2.30) 



Um 240 ^^ + "™ 10080 ' 



To calculate the corresponding Fade approximants for 



the second-order derivative, Um -, we begin with 



/.(-) 



B{E) 
HE) 



0{Ax^) 



(2.33) 



where J^{E) is given again by (2.21), and require that 
the approximation is fourth-oidei accurate or better. It 
follows that the operator B{E) must be symmetric in E, 
e.g. 



S(^). 



1 



aAx'^ 



{E 



2 , ^-2^ 



-p{E + E-^)+j 



and solve for 

HE)uli:^=B{E)um + 0{Ax''). 
We have 



(2.34) 
(2.35) 



B{E) Um = ^^ [2Um + 4u^ri:^Ax^ 

-^p[2um^uli:^Ax'^ul^^^^^- 



which gives the system of equations 

2/3 + 7= -2, 
a-/3 = 4, 
3a - /3 = 16 , 

with the solution 

a = 6, /3 = 2, 7 = 
Hence, we find 

,Ax2 
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and we can write 



,iii) 



Urrri^ U, 
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B2{E) 

_.Aviii)f^ 
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V840 rJ 24 



0(Ax«) , 



with 



B2{E) = 



6Ax2 



E' 



2E-^^2E-^ ^E- 



. (2.43) 



D. Fourth-order derivatives 



Because we chose to begin our derivation by improv- 
ing the third-order derivative, Um\ and both the finite- 
differences approximation for Um and Um already in- 
volve the entire subset, {x^, Xm±i,Xm±2}, it follows that 
we are limited to a second-order accurate approximation 
for the fourth-order derivative, Um . The error corre- 
sponding to the Pade approximant. 






0{Ax^) , (2.44) 



is obtained from the equation 

J^{E) ii^-) = Vr{E) u^ + (9(Ax2) . (2.45) 

Using F{E) from Eq. (2.21) and 



A ^2 A 4 



(2.46) 



find 

vt^ = ^M urn^u^::^ ^ + 0{Ax') . (2.47) 



HE) 



12 



III. APPROXIMATION SCHEMES 

Based on the above considerations regarding 
Pade approximants on the subset of grid points, 
{xm,Xm±i^Xm±2}^ it follows that wc cau always ob- 
tain a scheme that provides fourth-order accurate 
approximations for the derivatives Umj Um\ and Um • 
It is however possible to obtain approximants that 
mix fourth-order accurate approximations for two of 
these derivatives with a sixth-order accurate Pade 
approximant for the third one. We will discuss these 
special cases next, together with what may represent the 
"optimal" fourth-order approximation scheme. 

(6,4,4) scheme: This approximation scheme is an 
extension of the scheme introduced by Sanz-Serna et 
al. [29, 30] using a fourth-order Petrov-Galerkin finite- 
element method, and corresponds to choosing r = 30 in 
Eqs. (2.31) and (2.22). Then, we have 



120 , 6 = 26 , 



66, 



(3.1) 



and the coefficient of Ax^ vanishes in Eq. (2.31). There- 
fore, we obtain a sixth-order accurate approximation for 
the first-order derivative. 



/(O 



A2{E) 

-^[644] {E) 



(vii) _ 



5040 



• O(Ax^) , (3.2) 



a fourth-order accurate approximation for the second- 
order derivative. 



Xii) 



B2{E) 

-^[644] {E) 



(vi)'_ 



720 



• 0{Ax^) , (3.3) 



and a fourth-order accurate approximation for the third- 
order derivative, 

v^"^ = ^^^ u^+ut"^ ^ + 0( Ax«) , (3.4) 

>[644](^j 240 

where we introduced the notation 



-^[644] {E) 



120 



E'^ + 26£; -^m^ 26E-^ + E''^ 



(3.5) 



(4,6,4) scheme: The coefficient of Ax^ in Eq. (2.42) 
does not vanish for an integer value of r. To obtain a 
sixth-order accurate approximation for Um \ we require 
r = 180/7. Then, we have 



720 



a = 



b = 



152 



402 



c = 



(3.6) 



and we obtain a fourth-order accurate approximation of 
the first-order derivative, 



,(«) 



^'^^yu^+u(Z^^ + 0{Ax'), (3.7) 
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720 



a sixth-order accurate approximation of the second-order 
derivative, 



Xii) 



^'^^^ ■^i^+^i(^-)-ii-Ax^ + 0(Ax^), 



-^[464] {E) 



60480 



(3.^ 



and a fourth-order accurate approximation of the third- 
order derivative, 

v^"^ = ^^^ u^+ut-^ ^ + 0( Ax«) , (3.9) 

'>^[464](^j loO 

where we introduced the notation 



-^[464] {E) 



720 



lE'^ + 152£; + 402 + 152£;-^ + 7^"^ 

(3.10) 



(4,4,6) scheme: For r = 60, the coefficient of Ax^ van- 
ishes in Eq. (2.22) and we obtain a sixth-order accurate 
approximation for Um • We have 

a = 240 , 6 = 56 , c = 126 , (3.11) 

and we obtain a fourth-order accurate approximation of 
the first-order derivative. 
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FIG. 1. (Color online) Study of the K(2,2) (Rosenau-Hyman) compacton stability: For each numerical scheme we illustrate 
results at time t=175, the compacton propagation in its comoving frame with At=0.002 and Ax=0.1, 0.05, and 0.02. In all 
cases the radiation depicted here is a numerical artifact that is suppressed by reducing the grid spacing, Ax. This indicates 
that indeed the compacton is a stable solution of the K(2,2) model. 



a fourth-order accurate approximation of the second- 
order derivative, 
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and a sixth-order accurate approximation of the third- 
order derivative, 
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where we have introduced the notation 
1 



-^[446] (£^) = 



240 



£;^ + 56^ + 126 + 56^"^ +£;"^ . 

(3.15) 



This scheme is an extension of the scheme introduced 
first by Rus and Villatoro [33, 34]. 



(4,4,4) scheme: Finally, for the smallest value of r 
leading to integer positive values of a, 6, and c (i.e. r = 
5), we obtain 



a = 20, 6=1, c=16. 



(3.16) 



This gives a fourth-order accurate approximation of the 
first-order derivative, 



,(i) 



■^'^^^ u^+ui^^^+OiAx'), (3.17) 
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a fourth-order accurate approximation of the second- 
order derivative, 

n^ = ^^ u^-u(-^ §- Ax' + 0{Ax') , (3.18) 

and a fourth-order accurate approximation of the third- 
order derivative. 
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1 r. 
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IV. RESULTS 



. (3.20) 



To compare the quahty of the approximations dis- 
cussed above, we speciahze to the case of the K{p^p) 
equation. In a frame of reference moving with velocity 
Co, the K{p^p) equation reads 



du 



Co 



du 
dx 



duP 
dx 



d^uP 
dx^ 



= 0, l<p<3. (4.1) 



For p restricted to the interval 1 < p < 3, the K{p^p) 
equation allows for a compacton solution, with the simple 
form [32, 33, 40] 



Uc{x,t) 



a ' cos 



27 



(S^x.t) 



|e(x,t)|<^/(2/3), 

(4.2) 
where c is the compacton velocity and xq is the position 
of its maximum at t = 0, and we have introduced the 
notations ^(x, t) = x — xq — {c — co)t, and 



a 



2cp ^^_p 



1 



p+l 



2p 



7 



1 



p-1 



(4.3) 



Numerically, the lack of smoothness at the edge of 
the compacton introduces numerical high-frequency dis- 
persive errors into the calculation, which can destroy 
the accuracy of the simulation unless they are explicitly 
damped (see e.g. discussion in Ref. [25]). As such, we 
solve Eq. (4.1) in the presence of an artificial dissipation 
(hyperviscosity) term based on fourth spatial derivative, 
/i d^u/dx^^ and we choose ji as small as possible to reduce 
these numerical artifacts while not significantly chang- 
ing the solution to the compacton problem. We note 
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FIG. 2. (Color online) Comparison of compacton-stability 
results as a function of numerical scheme. Results are shown 
at t=175, for compacton propagation in its comoving frame 
with At=0.002 and Ax=0.1. 



nonetheless, that the addition of artificial dissipation re- 
sults in the appearance of tails and compacton amplitude 
loss. 

Let us consider now the numerical solution of Eq. (4.1) 
by means of the fourth-order accurate Fade approximants 
discussed here. In general, we can discretize Eq. (4.1) in 
space as 

"^"^ \coA{E) - ^iV{E)\um 



HE) 



At 



A{E) + C{E) {u. 



0, 



(4.4) 



We consider a uniform grid in the interval x G [0, L] 
by introducing the grid points Xm = mAx, with m = 
0, 1, • • • , M and the grid spacing Ax = L/M. In Eq (4.4), 
we assume that Um{t) obeys periodic boundary condi- 
tions, Uuif) = Uo{t). 

Following RV [ ] , we have numerically discretized the 
time-dependent part of Eq (4.4) by implementing both 
the implicit trapezoidal (Euler) and the implicit midpoint 
rule in time. Correspondingly, we need to solve the fol- 
lowing two approximate equations for Eq. (4.4): 



HE) 



,n+l 



At 

coAiE) - nV{E) 



(4.5) 



,n+l 



A{E) + C{E) 



«+i)^ + «)P 



corresponding to the trapezoidal rule, and 



HE) 



,n+l 



At 



(4.6) 



coA{E) - iiV{E) 
A{E) + C{E)'\ (- 



,,n+l 



n+1 



2 



= 0, 



corresponding to the midpoint rule. Here we have intro- 
duced the notations, u]^ = Um{tn) and u]^^ = Um{tn + 
At). 

In the following, we further specialize to the case of 
the i^(2, 2) equation (p = 2), which allows for the exact 
compacton solution 



Uc[X,t) = —COS 



X — {c — Co)t 



(4.7) 



in the interval \x — {c — co)t\ < 27r, where c is the velocity 
of the compacton. We note that in our simulations per- 
taining the i^(2, 2) compacton problem, we did not find 
any numerically-significant differences between stepping 
out the solution using the trapezoidal and the midpoint 
rules. This is consistent with the observation made by 
RV in Refs. [33, 34]. Therefore, in the following we only 
present results obtained using the trapezoidal rule. Im- 
plementing both methods is however important for the 
purpose of future simulations of compactons exhibiting 
higher-order nonlinearities, e.g. in the case of the L(r, s) 
compactons. 




FIG. 3. (Color online) Collision of two compactons with 
ci — 1 and C2 — 2. The simulation is performed in the co- 
moving frame of reference of the first compacton, i.e. cq — ci, 
using the (6,4,4) scheme and a hyperviscosity, /j, = 10~^. The 
collision is shown to be inelastic, despite the fact that the 
compactons maintain their coherent shapes after the collision. 




FIG. 4. (Color online) The first compacton (ci — 1) is "at rest" 
before the collision depicted in Fig. 3. As shown in the inset 
(with the same axis labels), after the collision the centroid of 
this compacton changes position and the compacton moves 
slowly consistent with a small change in amplitude due to 
hyperviscosity. 



A. Study of compacton stability 

To illustrate a numerical study of a compacton stabil- 
ity problem, we apply the Pade approximations discussed 
above to the case of the i^(2, 2) compacton defined in 
Eq. (4.7). The numerical compactons propagate with the 
emission of forward and backward propagating radiation. 
In Fig. 1, we illustrate results for each numerical scheme 
by depicting the numerically-induced radiation in the co- 
moving frame of the compacton (cq = c). Here we chose 
a snapshot at t=175 after propagating the compacton in 
the absence of hyperviscosity (/i=0) with a time step, 
At=0.002, and grid spacings, Ax=0.1, 0.05, and 0.02. 
We note that the amplitude of the radiation train is at 
least 7 orders of magnitude smaller than the amplitude of 
the compacton. Using the grid refining technique, we can 
show that indeed the radiation is a numerically-induced 
phenomenon. The noise is suppressed by reducing the 
grid spacing. Ax, indicating that the compacton (4.7) is 
a stable solution of the i^(2, 2) equation. 

For any numerical study of compacton stability and 
dynamical properties, it is important to reduce as much 
as possible the numerically-induced radiation. It is de- 
sirable to minimize three characteristics of the radiation 
train: (i) the length of the radiation train, in order to 
avoid a wrap around of the solution as a result of the 
periodic boundary conditions constraint, (ii) the ampli- 
tude of the radiation train, which should be minimized 
in order to better differentiate numerical artifacts from 
physics, and (iii) the amplitude at the leading edge of 
the radiation train, which seems related to the suscep- 
tibility of the numerical approximation to instabilities 
arising particularly in dynamical studies. Large ampli- 
tudes at the leading edge of the radiation train lead to 
the need for large values of the hyperviscosity parameter, 
/i, in order to overcome these instabilities. 

To study the quality of our Pade approximations, in 
Fig. 2 we compare the radiation results at t=175 fol- 



T^^ 0.0 



y: 0.0 



'^- 0.0 



X 0-0 



;^ 0.0 



^: 0.0 



><^ 0.0 



y, 0.0 



Z: 0.0 



X 0.0 



;^ 0.0 





75 80 85 90 95 100 105 110 115 120 



FIG. 5. (Color online) Dynamics of the zero-mass "ripple" 
with shock created as a result of the collision depicted in 
Fig. 3. 



FIG. 6. (Color online) Comparison of the Fade numerical 
schemes in the context of the ripple created as a result of 
the collision depicted in Fig. 3. In the upper panel we illus- 
trate the ripple calculated for t=80 using the (6,4,4) scheme, 
whereas in the bottom panel we illustrate the differences be- 
tween results obtained using the other schemes and the (6,4,4) 
scheme. All simulations were performed using a hyperviscos- 
ity, /i = 10-^ 



and third-order spatial derivatives. It appears, that at 
least for the i^(2, 2) compacton, an improved first-order 
derivative approximation, i.e. the (6,4,4) scheme, leads to 
a shorter radiation train than in the case of the (4,4,6) 
scheme, which improves the quality of the third-order 
derivative [ ]. However, the amplitude of the radiation 
train wave in the (4,4,6) scheme is comparable with the 
train amplitude in the (4,4,4) scheme, and smaller than 
the train amplitude in the (6,4,4) scheme, which indi- 
cates that it is important to improve the numerical ap- 
proximation of the third-order spatial derivative in the 
i^(2,2) equation. The i^(2,2) equation does not fea- 
ture a second-order spatial derivative, and the (4,6,4) 
scheme behaves as a tradeoff between the (6,4,4) and 
(4,4,6) schemes: the radiation train is shorter in the 
(4,6,4) scheme, but the amplitude of the train is com- 
parable with that in the (6,4,4) scheme and larger than 
in the (4,4,6) scheme. Finally, with respect to the ampli- 
tude at the leading edge of the radiation train, the (4,4,6) 
scheme is the best and the (4,4,4) scheme is the worst and 
likely will require the largest hyperviscosity parameters 
in dynamical problems. 



B. Study of compacton dynamics 



lowing the compacton propagation with At=0.002 and 
Ax=0.1 in the comoving frame of the compacton. We 
notice that the spatial extent of the radiation train is 
minimum in the optimal (4,4,4) scheme. Furthermore, 
the K{p^p) equation, Eq. (4.1), depends only on first- 



It is generally accepted that the RH compactons have 
the soliton property that after colliding with other com- 
pactons they reemerge with the same coherent shape. 
However, unlike in soliton collisions, the point where 
the compactons collide is marked by the creation of low- 
amplitude compacton-anticompacton pairs [ ]. De Fru- 
tos et at. showed [ ], and RV confirmed recently [33], 
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FIG. 7. (Color online) Similar to Fig. 6. Here we compare 
the numerical schemes in the context of the shock compo- 
nents observed in the ripple created as a result of the collision 
depicted in Fig. 3. 



tions were performed using a hyperviscosity, /i = 10~^. 
Similarly, in Fig. 7 we compare our numerical schemes in 
the context of the shock components observed when the 
ripple switches from negative to positive values. 

Based on the results depicted in Figs. 6 and 7, we 
observe that, independent of the numerical scheme, the 
largest errors occur at the end of the ripple in the direc- 
tion of its propagation (see Fig. 6), and the errors are 
very similar in the area of the two shock components 
(see Fig. 7). Furthemore we note that the results of 
schemes (6,4,4) and (4,4,4) are very similar. Given that 
based on our stability studies we concluded that scheme 
(4,4,4) provides the most accurate set of Fade approxi- 
mants for the i^(2, 2) problem, this leads us to use the 
(6,4,4) scheme as the reference for this comparison. Fi- 
nally, the (6,4,4) results are closer to the (4,4,6) results 
than they are to the (4,6,4) results, which seems to indi- 
cate that for dynamical i^(2, 2) problems (4,6,4) approx- 
imation scheme fares the worst, as the i^(2, 2) equation 
does not depend on second-order spatial derivatives. 



that shocks are generated during compacton collisions. 
Shocks are also generated when arbitrary initial "blobs" 
decompose into a series of compactons [ ] . These shocks 
offer an ideal setting to compare numerical approxima- 
tions such as the Fade approximants discussed here. 



1. Pairwise interaction of compactons 

In this scenario, we consider the collision between two 
compactons (4.7) with velocities Ci = 1 and C2 = 2. In 
Fig. 3, we depict a series of snapshots of this collision 
process. The compactons are propagated in the comov- 
ing frame of reference of the first compacton, i.e. cq = ci, 
using the (6,4,4) scheme and a hyperviscosity, fi = 10~^. 
The collision is shown to be inelastic, despite the fact 
that the compactons maintain their coherent shapes af- 
ter the collision. The first compacton is "at rest" before 
the collision occurs. After the collision, this compacton 
emerges with the centroid located at a new spatial loca- 
tion, as illustrated in Fig. 4 . The inset in Fig. 4 shows 
the compacton moving slowly after collision, consistent 
with a small change in its amplitude. 

The collision process depicted in Fig. 3 gives rise to a 
zero-mass ripple with a shock when the "ripple" switches 
from negative to positive values (see Fig. 5). A small 
change in the shock amplitude is noticed and is due to 
the presence of hyperviscosity. These shock components 
were first noted in Ref. [ ] , and were shown to be robust 
with respect to the numerical approximation in Ref. [34] . 
In Fig. 6 we compare results obtained with our four Fade 
approximants schemes in the context of the ripple created 
as a result of the collision depicted in Fig. 3. In the up- 
per panel of Fig. 6 we show the result obtained using the 
(6,4,4) scheme at t=80, whereas in the bottom panel we 
illustrate the differences between results obtained using 
the other schemes and the (6,4,4) scheme. All simula- 



2. Dynamics with arbitrary initial conditions 



Following RV 
"blob" given 



we consider the time evolution of a 



u{x^O) 



y COS^ ^ 



3 ' 

^ cos^ — 



^^° for 150 - 27r < X < 150 . 
for 150 < 160 , 



160 



for 160- < X < 160 



•27r. 
(4.8) 

In Fig. 8 we illustrate the dynamics of this blob de- 
composition, as calculated using the (6,4,4) scheme and 
a hyperviscosity, /i = 10~^, and show that the blob 
evolves into two compactons and a ripple featuring a set 
of compacton-anticompacton pairs. Similar to the colli- 
sion problem, the ripple has positive- and negative- value 
components separated by a shock. 

As surmised following our compacton stability study, 
the radiation train corresponding to the (4,4,4) scheme 
has a higher amplitude at the leading edge, which makes 
this scheme more susceptible to instabilities than the 
other three schemes. This undesirable feature of the 
(4,4,4) scheme is illustrated in the case of this blob de- 
composition. In Fig. 9, we depict the last three steps in 
the simulation and compare results obtained using the 
(6,4,4) and (4,4,4) schemes. These results show how the 
(4,4,4) simulation of the "blob" decomposition becomes 
unstable and crashes corresponding to a hypeviscosity, 
jii = 10~^, whereas the (6,4,4) simulation does not. The 
(4,4,4) simulation becomes stable if we increase the hy- 
perviscosity (see Fig. 10), but the instability cannot be 
removed by reducing the grid spacing (not shown). 

With the exception of the instability developed in the 
(4,4,4) scheme, the comparison of the four Fade schemes 
reveals a situation very similar to the case of the pairwise 
compacton collision and results are illustrated in Fig. 11 
at t=15.5 [just prior to scheme (4,4,4) becoming unstable] 
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FIG. 8. (Color online) Dynamics of a "blob" decomposi- 
tion into two compactons and a ripple featuring a set of 
compacton-anticompacton pairs. Similar to the collision 
problem, the ripple has positive- and negative- value com- 
ponents, separated by a shock front. The simulation was 
performed using the (6,4,4) scheme and a hyperviscosity, 
/i=10-^ 



and Fig. 12 at t=60. By comparing the four approxima- 
tion schemes in the context of the shock front formed 
when the ripple switches from positive to negative val- 
ues, we find that results obtained using schemes (4,4,4) 
and (6,4,4) are very close to each other and the (6,4,4) 
results are closer to the (4,4,6) results than they are to 




FIG. 9. (Color online) Autopsy of an instability: The simula- 
tion of the "blob" decomposition depicted in Fig. 8, becomes 
unstable and crashes when performed using the (4,4,4) scheme 
with a hyperviscosity, /i = 10~^. Here we illustrate the last 
three steps in the simulation by comparing the results ob- 
tained using the (6,4,4) and (4,4,4) schemes. The insets show 
a magnified region around x = 150 for clarity. 



the (4,6,4) results. Unfortunately, the (4,4,4) scheme re- 
quires a larger hyperviscosity to smooth out the noise 
at the shock front. Our study suggests that the (6,4,4) 
scheme is the best alternative to scheme (4,4,4). 



V. CONCLUSIONS 

To summarize, in this paper we presented a system- 
atic approach to calculating higher-order derivatives of 
smooth functions on a uniform grid using Fade ap- 
proximants. We illustrated this approach by deriving 
higher-order approximations using traditional second- 
order finite-differences formulas as our starting point. 
We proposed four special cases of the fourth-order Fade 
approximants and employed these schemes to study the 
stability and dynamical properties of K{p^p) Rosenau- 
Hyman compactons. This study was designed to estab- 
lish the baseline for future studies of the stability and 
dynamical properties of L{r^s) compactons, Eq. (1.2), 
that unlike the RH compactons are derivable from a La- 
grangian ansatz. The L(r, s) compactons feature higher- 
order nonlinearities and terms with mixed-derivatives 
that are not present in the K{p,p) equations, and hence 
the numerical analysis of their properties is expected to 
be considerably more difficult. Finally, we intend to ap- 
ply these schemes to the study of PT-symmetric com- 
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FIG. 10. (Color online) Smoothing effect of the hyperviscos- 
ity: The simulation of the "blob" decomposition performed 
using the (4,4,4) scheme becomes unstable for a hyperviscosity 
jj — 10~^. The instability cannot be removed by reducing the 
grid step (not shown) . Instead one must increase the value of 
the hyperviscosity, which results in a smoothing of the noise 
at the shock front. 




FIG. 11. (Color online) Comparison of the numerical schemes 
discussed in the context of the shock front formed when the 
ripple switches from negative to positive values. In the up- 
per panel we depict the result at t=15.5 obtained using the 
(6,4,4) scheme, whereas in the bottom panel we illustrate the 
differences between results obtained using the other schemes 
and the (6,4,4) scheme. All simulations were performed using 
a hyperviscosity, /x = 10~^. Shortly after t — 15.5 the simu- 
lation performed using the (4,4,4) scheme becomes unstable. 



pactons [2>^]. 

Based on our compacton stability study, we conclude 
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FIG. 12. (Color online) Similar to Fig. 11. Here, results 
are presented for t=80. The simulation performed using the 
(4,4,4) scheme became unstable shortly after i — 15.5 and is 
not represented here. 



that none of our four Fade approximations appears as a 
clear winner in the context of the three minimization cri- 
teria for an optimal numerical discretization of the com- 
pacton problem: length and amplitude of the radiation 
train, and amplitude at the leading edge of the radia- 
tion train. The (4,4,4) scheme features the shortest and 
smallest amplitude of the radiation train, but exhibits the 
largest amplitude at the leading edge of the train. The 
(4,4,6) scheme has the smallest amplitude at the leading 
edge of the train and the amplitude of the train is com- 
parable otherwise with that in the (4,4,4) scheme, but 
the length of the radiation train in the (4,4,6) scheme is 
the largest of the four approximations considered here. 
Hence, the (4,4,6) scheme will probably be most useful 
in studies of the dynamical properties of the system and 
deal best with shock-type problems, but will require the 
largest model space (largest value of V) to avoid the un- 
desired wrap around effect of the solution due to the pe- 
riodic boundary conditions constraint, which in turn will 
make grid refining studies difficult in the (4,4,6) scheme, 
because of physical constraints in allocatable memory 
and CPU wall time. 

The best compromise may be provided by the (6,4,4) 
approximation scheme. Our dynamical simulations indi- 
cate that the results obtained using this method closely 
resemble those obtained using the (4,4,4) scheme, and the 
(6,4,4) scheme is less susceptible to instabilities, at least 
in the particular cases discussed here, than the (4,4,4) 
scheme. 
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Appendix A 

In this appendix we discuss the derivation of fourth- 
order accurate approximations for the first- and second- 
order derivatives of a smooth function. 

To derive a fourth-order accurate approximation of the 
first-order derivative, we begin by introducing the oper- 
ator 
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and ask that the following relation is fulfilled: 
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To satisfy the requirement of a fourth-order accurate ap- 
proximation for the first-order derivative i^^ , we solve 
the system of equations 



a - 2/3 = 4 , 
P= -8, 

and obtain the solution 

a = -12, /3 = -8. 



(A4) 
(A5) 



(A6) 
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This gives 



We have 
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with 



2u^ + u^^^^Ax'+ui^^^ + ^ 



which gives the system of equations 



-7^7 



.} 



(All) 



Ai{E) = 



1 



12Aa; 



E^ -8E + 8E-^ - E- 



■ (A8) 



2/3 + j= -2, 


(A12) 


a- /3 = 4, 


(A13) 


3a- 13 = 0, 


(A14) 



Similarly, to derive a fourth-order accurate approxi- ^ith the solution 
mation for the second-order derivative we introduce the 
operator 



Bi{E)u, 



aAx'^ 



[E^^E-^)^p[E^E-^)^^ 



Urn • 

(A9) 



a =-2, /3=-6, 7=10. 
Hence, we obtain 
1 



(A15) 



Bi{E) 



2Aa;2 



E^ -6E+10- 6E-^ + E-^ 



(A16) 



and seek a, /3, and 7 such that 

ul^'^=Bi{E)u^ + 0{Ax^). 



and 



,ri\ 



(AlO) 



4:^ = Bi (E) um + u^') - Ax^ + 0(Ax6) . (A17) 



